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SYMBOLS 


A mass-flux constant, equation (1 ) 

a sound speed at fixed composition, cm/sec 

B momentum-flux constant, equation (2) 

C total enthalpy constant, equation (3) 

Cpj? constant-pressure specific heat at frozen composition, cal/initial mole-°K 
D dissociation energy (1 17,960 cal/mole for 0 2 ); also, total derivative 

G nonequilibrium characteristics term, equation (25), g/cm 3 -sec 

h static enthalpy per initial mole of molecules, cal/mole 

K c equilibrium constant, moles/cm 3 

krf dissociation rate, cm 3 /mole-sec 

k r recombination rate, cm 6 /mole 2 -sec 

O oxygen atom 

0 2 oxygen molecule 

P pressure, atm, or dynes/cm 2 

p partial pressure, atm 

R universal gas constant, 1 .987 cal/mole-°K, or 82.06 cm 3 -atm/mole-°K 

R r specific gas constant, 2.598X10 6 ergs/g-°K for0 2 

S entropy per unit mass, cal/g-°K 

T temperature , ° K 

t time, sec 

U velocity of incident shock wave relative to end wall, cm/sec 

u velocity of gas behind reflected shock wave relative to end wall, cm/sec 

V velocity of reflected shock wave relative to end wall, cm/sec 

iii 
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v velocity of gas behind incident shock wave relative to incident shock, cm/sec 

W molecular weight of diatomic gas (32 g/mole for 0 2 ) 

w mass rate of production of a chemical species, g/cm 3 -sec 

x distance from shock-tube end wall, cm 

a degree of dissociation (equals, p Q /p) 

P arbitrary multiplier of recombination rate 

5 derivative along characteristic line 

Orf characteristic temperature of dissociation (58,960° K for 0 2 ) 

8 V characteristic temperature of vibration (2,274° K for 0 2 ) 

I u chemical potential, cal/g 

p mass density, g/cm 3 

r characteristic time 

[ ] concentration, moles/cm 3 

Subscripts 

eq equilibrium 

k specific chemical species 

1 ahead of incident shock wave 

2 behind incident shock wave 

3 behind reflected shock wave 

3 o behind reflected shock wave at instant of reflection 

2 1 general spatial point behind incident shock wave 



NUMERICAL SOLUTIONS OF SEVERAL REFLECTED SHOCK-WAVE 


FLOW FIELDS WITH NONEQUILIBRIUM 
CHEMICAL REACTIONS 

Ronald K. Hanson,* Leroy L. Presley, and Elsie V. Williams 
Ames Research Center 


SUMMARY 


The method of characteristics for a chemically reacting gas is used in the construction of the 
time-dependent, one-dimensional flow field resulting from the normal reflection of an incident 
shock wave at the end wall of a shock tube. Nonequilibrium chemical reactions are allowed behind 
both the incident and reflected shock waves. All the solutions are evaluated for oxygen, but the 
results are generally representative of any inviscid, nonconducting, and nonradiating diatomic gas. 
The solutions clearly show that (1) both the incident- and reflected-shock chemical-relaxation 
times are important in governing the time to attain steady-state thermodynamic properties; and 
(2) adjacent to the end wall, an excess-entropy layer develops wherein the steady-state values of 
all the thermodynamic variables except pressure differ significantly from their corresponding 
Rankine-Hugoniot equilibrium values. 


INTRODUCTION 


The use of shock tubes for studies of high-temperature gasdynamics and nonequilibrium 
chemical kinetics is well established. Many of these studies consisted of observations of the gas 
heated by the reflection of an incident or primary shock wave at the end wall of a shock tube 
(Johnson and Britton, ref. 1). Reflected-shock testing has two major advantages: First, the 
temperature behind the reflected shock wave is increased to roughly twice that behind the incident 
shock wave. Second, the gas behind the reflected shock wave is nearly stationary in laboratory 
coordinates; thus greater time resolution can be achieved for studying nonequilibrium processes, for 
example, than is obtainable behind shock waves at the same thermodynamic conditions. 

However, several phenomena occur behind a reflected shock wave that can limit its usefulness 
in producing a well-defined sample of hot test gas. For the most part, these phenomena give rise 
to uncertainties in the thermodynamic state of the test gas. Three phenomena occur as a consequence 
of the side-wall viscous boundary layer that develops behind the incident shock wave. First, the 
presence of the boundary layer induces spatially nonuniform properties within the core of the 
flow field behind the incident shock wave, thus causing time-dependent upstream properties for 
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the reflected shock wave. From their analyses of this problem, Kamimoto et al. (ref. 2) and 
Rudinger (ref. 3) have concluded that the nonuniform flow is responsible for the slow pressure 
rise often observed on the end wall of a shock tube after shock reflection. The second phenomenon 
is the initiation of a complicated shock-wave-boundary-layer interaction as the reflected shock 
propagates back through the side-wall boundary layer. This interaction is more intense in gases of 
greater molecular complexity. Mark (ref. 4) and Strehlow and Cohen (ref. 5) have observed the 
bifurcation resulting from the interaction and concluded that the effect is negligible for monatomic 
gases and is sufficiently small for diatomic gases to permit kinetics studies immediately adjacent to 
the end wall. In an interferometric study of strong shocks in dissociating nitrogen, Hanson 
(unpublished results (1968)) observed increasing bifurcation with increasing shock speed and 
decreasing test pressure, and found that the flow near the end wall remained suitable for kinetics 
studies only for the short times prior to serious bifurcation effects. Third, the side-wall boundary 
layer causes the incident shock wave to be curved. On the basis of experimental data, Lin and Fyfe 
(ref. 6) argued that the center displacement of the incident shock wave is directly proportional to 
(r/pi ) 1,/2 , where r is the shock -tube radius and Pi is the initial density. The effects of a curved 
incident shock wave on a reflected-shock-wave flow field, however, are not well understood. 

A fourth phenomenon that complicates the reflected-shock-wave flow field arises from the 
thermal boundary layer that forms on the end wall of the shock tube. Sturtevant and Slachmuylders 
(ref. 7) have studied, both theoretically and experimentally, the effect of the end-wall thermal 
boundary layer on the trajectory of the reflected shock wave; and Goldsworthy (ref. 8) has obtained 
expressions for the flow-field perturbations resulting from the boundary-layer formation. Generally, 
these results show that if the initial tube pressure is greater than about 1 torr, the effects of the heat 
transfer on the flow field are not significant at distances greater than 1 mm from the end wall, 
provided observation times are less than a few hundred microseconds. 

Since all the aforementioned phenomena result from viscous or heat-conduction effects, the 
controlling factors, if not the phenomena themselves, are well understood. These complications 
can thus be minimized by appropriate selection of experimental technique, test conditions, and 
observation interval. However, an additional problem exists for studies involving nonequilibrium 
chemical reactions because an adequate theory for predicting the entire, nonequilibrium flow field 
in even an inviscid and nonconducting gas has not been available. This latter aspect of shock-wave 
reflection is the subject to which this paper is addressed. A simple model for shock-wave reflection 
in a relaxing gas has been presented by Hanson (ref. 9), but results based on that model apply only to 
the layer of gas immediately adjacent to the end wall. The problem here is to determine the entire 
reflected-shock flow field when nonequilibrium chemical reactions are present behind both the 
incident and reflected shock waves. The method of characteristics for a one-dimensional reacting 
gas is used in obtaining numerical solutions. 

The closely related problem of reflected-shock flow fields with vibrational nonequilibrium 
was first discussed qualitatively by Baganoff (ref. 10) with particular emphasis on the vibrational 
relaxation of carbon dioxide. More recently, Johannesen et al. (ref. 1 1) have obtained quantitative 
results for the vibrational relaxation of carbon dioxide behind reflected shock waves by means of the 
method of characteristics. Many of the general features of these works also appear in the present 
results, although the emphasis of this work is on defining quantitatively the effect of nonequilibrium 
chemical reactions on a reflected-shock flow field. A simple chemical model has been used to 
show clearly the coupling between the fluid mechanics and chemistry. Once the basic features of a 
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nonequilibrium reflected-shock flow field are understood, a more complicated chemical model can 
readily be implemented. 


METHOD OF ANALYSIS 


The problem of finding a solution for a reflected-shock flow field with nonequilibrium chemical 
reactions introduces several phenomena that are not considered in the usual equilibrium solutions 
of reflected-shock properties. For example, nonequilibrium flow behind the incident shock wave 
causes curvature of the particle paths in the xt plane (see fig. 1 ), and results in a spatial gradient in 
all the fluid and thermodynamic properties behind the incident shock wave. The reflected shock 
wave must propagate back into this nonuniform region; hence the properties upstream of the 
reflected shock are continually changing. As a result of both the time-varying upstream properties 
and the nonequilibrium chemical reactions occurring behind the reflected shock wave, the velocity 
of the reflected shock wave is not constant (fig. 1). In addition, the usual boundary condition of 
zero flow velocity behind the reflected shock applied for an ideal gas now holds only at the end 
wall, and the particle paths behind the reflected shock wave are curved. These factors complicate 
the flow-field calculations by preventing selection of a coordinate system in which the flow across 
the reflected shock wave appears steady, as is possible with incident-shock-wave calculations. 

Since the relaxation region behind the incident shock wave is independent of the reflected 
shock wave, it is possible to separate the overall solution for shock reflection into two portions: 
(1) flow behind the incident shock wave, and (2) flow behind the reflected shock wave. Several 
assumptions regarding the mathematical model used to obtain solutions are: 

1 . The gas is a mixture of thermally perfect gases. 

2. The vibrational energy is in equilibrium at the local translational temperature and is given 
by a nontruncated harmonic-oscillator model. 

3. All molecules and atoms remain in their ground electronic states. 

4. The reaction under study is the dissociation of molecular oxygen, 

k d 

0 2 + A ^20 + X 

kr 

where X denotes an arbitrary collision partner that may be either molecular or atomic oxygen. 

5. The ratio of the dissociation and recombination rate constants is the equilibrium constant 
evaluated at the local translational temperature. 

6. The flow is one-dimensional in a constant-area tube. 

7. All transport processes are neglected. 

8. The end wall is chemically inert. 
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These assumptions greatly simplify the study without changing the principal features of the coupling 
between nonequilibrium chemical reactions and fluid mechanics in a reflected-shock flow field. 


Incident-Shock Solution 

The solution of incident-shock-wave flow fields with finite-rate chemistry is straightforward 
(ref. 12), primarily because the flow is steady in a shock-fixed coordinate system. The set of equa- 
tions describing the steady, one-dimensional flow of a dissociating diatomic gas behind a normal 
shock wave are as follows; in shock-fixed coordinates: 

Conservation equations. 

Mass: 


P\U = ( pv) 2 i = A 

Momentum: 


p ! + Piie = p 2i + ( P v 2 u = b 

Energy: 

h l + (1/2 )IP = h 2i + (l/2)v 2 j 2 = C 


( 1 ) 

( 2 ) 

(3) 


Equations of state. 
Thermal : 


P = (1 + a)~pT (4) 

or 

P = (1 +ci)R'pT (5) 


Caloric: 


Rate equation. 


h = 


7 + 3 a 


RT + 


(1 - a) R6 V 
e &v/T - i 


+ ocD 


D cx 


( 6 ) 


(7) 


where 


k r = $AT m 


( 8 ) 


4 



k-d ~ K c k r 


(9) 


K c = 0.666 T l/2 e 9d/T (\-e~ 6v/T ^ 
The total molar concentration [2f] is simply 



while the molecular molar concentration is 


[ 0 2 ] 


1 - a 
1 + a 


[X] 


and the atomic molar concentration is 


[O] 


2a 

1 + a 


IX) 


( 10 ) 

(ID 


( 12 ) 


(13) 


The quantities A, B, and C are constants for each flow case, readily evaluated for specified initial 
conditions and shock speed. The rate constant k r actually used in the numerical calculations was 
chosen subjectively on the basis of several experimental studies (see, e.g., Green and Toennies, 
ref. 13). The values of A and m selected are 7.29X1 0 2 1 (cm 6 /mole 2 -sec) and -2, respectively; 
|3 is a multiplier for conveniently varying the magnitude of the rate constant. The simplified 
expression for the equilibrium constant K c is given in appendix A. Equations (4) through (13) 
also apply to the flow behind a reflected shock wave. 


The above set of equations is solved at each point 2 i (i = 0, 1, 2, . . . ) in the incident-shock 
flow field, subject to the following boundary conditions: (1 ) no dissociation occurs across the shock 
front, that is, the chemistry is frozen 1 at a = 0; and (2) equilibrium is reached when Da/Dt goes 
to zero. The procedure makes use of two independent equations for the gas velocity v 2 i obtained 
from equations (1), (2), (3), and (5): 

"(t)[(x) 2 (14) 


and 


v 2 j = [2(C-/2 2i )] l/2 (15) 

Immediately behind the shock, at point 20, the value of a is zero, and T 20 may be found by 
varying the temperature until the same value of v 20 is obtained in equations (14) and (15). 
Subsequent flow-field points are fixed by selecting the value of T 2 i according to 

T 2i = T 20 - iAT (16) 


1 In this paper, “frozen” chemistry refers to fixed composition (fixed a). 
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where AT is a computer input quantity usually taken to be 100° K. The variable <x 2 i is found by 
simple iteration using equations (14) and (15) with the known value of T 2 \. It is clear that for any 
given shock wave, once a single thermodynamic quantity is known, at any flow-field point, all the 
remaining variables can be found directly. 

The position of each flow-field point is still unknown, however, and must be found by 
numerical integration of the rate equation. The increment of length between adjacent flow-field 
points is given by 


(A x) 2i = 


<*21 0^2 (/- l ) 

(da/dx ) av 


(17) 


where 



(18) 


and the values of da/dx are computed from equation (7). (Note: Da/DT = v 2 i(da/dx) 2 i for the 
steady flow behind the incident shock wave.) The total distance of point 2 i from the incident 
shock front is 


x 2i = ^ (Ax) 2 / (19) 

7=1 

In practice, calculations continue so long as da/dx is positive; when da/dx goes negative, 
equilibrium has been passed, and the program goes back at least two flow-field stations and 
computes forward with smaller temperature steps (generally 10° K)until equilibrium is again passed. 
The program then backs up again and computes forward with a final, yet smaller temperature-step 
size (generally 2° K). This small-step portion of the program is necessary to minimize errors in 
computed relaxation length near equilibrium. 


The actual equilibrium conditions can be established independently of the rate equation, of 
course, from equations (1) through (6) and the law of mass action, 


K c = 


4a 2 


1 - a 2 


[X] 


4a 2 p 
1 - a W 


( 20 ) 


For a homonuclear diatomic gas, the appropriate equations to solve are: 


A 2 



j \j-_a W 
| a 2 4 


K c 


R'T{\ +a) 


and 



( 21 ) 
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( C-h ) 


( 22 ) 



The equilibrium constant K c is a function of T alone (eq. (10)) and h is a function of T 
and a, so that simultaneous solution of equations (21) and (22) yields the equilibrium values of T 
and a. behind the incident shock wave. Other equilibrium property values follow from the con- 
servation and state relations. 

Flow-field properties throughout the incident-shock nonequilibrium zone are stored in 
computer memory since they provide the time-varying upstream conditions for the reflected shock 
wave. 


Reflected-Shock Solution 

The entire flow field behind the reflected shock wave is constructed using the method of 
characteristics. A derivation of the appropriate characteristic equations is given in appendix B 
(see also Ferri (ref. 14) and Chu (ref. 15)); only the final form of the equations is used here. In 
addition, since the complete calculation procedure is very lengthy, only a brief description of the 
logic of the numerical solution is presented here; a more detailed discussion is provided in appen- 
dixes C, D, and E. 

Basic equations- The equations of motion for a reacting gas, as for a nonreacting gas, can 
be simplified to ordinary differential form along the characteristic directions specified by dx/dt 
= u ± a. For flow in chemical nonequilibrium, of course, the appropriate sound speed a is that at fixed 
composition (see ref. 16, p. 300). The ordinary differential equations for the change in pressure and 
velocity along the characteristic Mach lines are . 

(1 l/p a 2 )(5 + P/5/ ) + (l/a)(8 + u/8t) = G/p (23a) 


and 


(l/pa 2 )(8~P/8t) - (l/a)(8~u/8t) = G/p (23b) 

where 5 + /5/ and 8 /8t correspond to differentials taken along right- and left-running charac- 
teristics, respectively; that is, 


5 + 8/ = 9/9/ + (w+fl) 9/9x (24a) 

and 

5~/5/ = 9/9/ + (u-a) b/dx (24b) 

Equations (23a) and (23b) form a convenient basis for numerical calculations and, with the 
differential expressions along the streamline characteristics, are known as the “compatibility 
relations.” Note that the compatibility relations given above are identical to those for a nonreacting 
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gas except that the term G/p is added to account for the effects of chemical nonequilibrium. The 
term G is given by (appendix B) 


3 D__ V 1 

2 r T(e e ” /T -l) 1+ “J 

where the specific heat at constant pressure and composition, is defined by 

/dh\ = / 7 + 3a\ (1 ~ a)R6 „ 2 / v/r 

C pf \dT / P'Ot \ 2 /* + T2^vlT_ x y 



and wo , the net mass rate of production of atoms, is given by 

w 0 = Wk r (K c [ 0 2 ][X) - fO] 2 m) 
The sound speed at fixed composition, a , is given by (ref. 16) 

(1 +a)R'T 


,2 = 


/a 


1 


S,a 

C pf-R(\ +a)j 


(25) 


(26) 


(27) 


(28) 


Consideration of nonequilibrium reacting flows complicates the method of characteristics in 
that changes in concentration and entropy occur along particle paths. A determination of the 
conditions at a given point in a flow field must account for these variations as well as the changes in 
pressure and velocity prescribed by equations (23). Along a particle path (the same as a streamline 
in the xt plane), differences in composition and entropy are given by 


8a = 

St p 


(29) 


and 


f < 30 > 

k 

where k represents a specific chemical species, in this case O or 0 2 . 

The quantity wq is given by equation (27), and atom conservation requires 

w 0 = -wq 2 (31) 
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The chemical potentials Mfc appearing in equation (30) are given by (appendix A) 



(32) 

(33) 


where po and po 2 are the partial pressures of oxygen atoms and molecules, respectively. 

Logic of numerical solution— The present solution, like all characteristics solutions, is 
initiated from a set of points along a curve where flow-field properties are known, such as points 1 , 
2, and 3 of figure 2. First the location and then the flow-field properties of new points are computed 
along subsequent characteristics. In the present work, points are computed successively along 
left-running characteristics, as shown, so that point 4 is computed following point 3, point 5 
after point 4, and so on. The location of a new point follows from the intersection of characteristic 
rays from two previous points (exceptions are end-wall points and shock-wave points, which are 
found with only one characteristic ray). The properties at each new point are computed in two 
steps. First, the pressure and particle velocity are found by an appropriate simultaneous solution of 
equations (23a) and (23b); only one of these equations together with the proper boundary 
conditions is required for shock-wave and end-wall points. Second, the remaining properties at the 
new point are found by integrating equations (29) and (30) along the particle path that passes 
through the new point and originates from an interpolated point a between two known points on 
the preceding characteristic (fig. 2). The solution continues by successive computation of rays. 
Details of the network construction technique are discussed in appendix C. 

The conditions along an initial ray — for example, points 1, 2, and 3 — are found from a 
starting approximation. First, the conditions at the point of reflection, point 0, are calculated 
assuming chemistry frozen at a = 0 across the reflected shock wave. Next, the properties at 
point 1 are established by assuming that pressure and enthalpy remain constant for a very short 
time along the end wall. Finally, the conditions at the boundary, point 3, are found from a simul- 
taneous solution of equation (23b) and the appropriate shock-jump relations. Conditions at point 
2 are obtained by simple linear interpolation between the conditions at points 1 and 3. A more 
complete description of the network-starting technique is given in appendix D. 

The brief discussion here pertains to a fixed number of points along each ray. Additional 
points are added (by linear interpolation) along each successive ray to keep the mesh sufficiently 
small to assure accuracy. Furthermore, to obtain good accuracy and stability in the region immedi- 
ately behind the reflected shock wave, where very large property gradients are present, a fine mesh 
is incorporated. The number of points in this fine mesh is variable and depends on the magnitude 
of the property gradients. Further details of the network-generating scheme are provided in 
appendix E. To further minimize errors caused by large property gradients, the solution at each 
point in the entire flow field is found by use of an iterative procedure. The properties first calcu- 
lated for a new point are averaged with those of the known or reference points, the average values 
are assigned to the coordinates of the reference points, and the new point is recomputed. In addition, 
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equations (29) and (30) are integrated along each particle path by assuming a linear variation in 
pressure. Each major interval, such as from point a to point 5 of figure 2, is divided into 20 small 
intervals for the numerical integration (appendix C). 


RESULTS AND DISCUSSION 


Five solutions of reflected-shock flow fields have been obtained; the initial conditions for 
these solutions are summarized in table 1 . Detailed results of case 1 are discussed first to establish 

the general character of reflected-shock 


TABLE 1.- SUMMARY OF COMPUTER SOLUTIONS; 
= 294° K, ALL CASES 


Case 

u, 

km/sec 

P i, 

atm 

P 

1 

3.05 

2.63X10 -3 

1 

2 

3.05 

13.15X1 0 — 3 

1 

3 

3.05 

2.63X10 -3 

5 

4 

3.96 

2.63X10 -3 

1 

5 

2.01 

2.63X10 -3 

1 


flow fields with nonequilibrium chemical 
reactions. Results from cases 2 through 5 
are presented to show the effect of changing 
the initial pressure (case 2), the reaction- 
rate constant (case 3), and the incident- 
shock velocity (cases 4 and 5). Preliminary 
results are given in reference 17. 

The results of the computer program 
are displayed in three different coordinate 
planes: incident-shock flow-field profiles 
as a function of distance behind the incident 
shock; end-wall properties (x = 0) as a 
function of time; and flow-field properties 
as a function of x for discrete times. 


Case 1 

Profiles of flow-field properties for the incident-shock solution of case 1 are shown in 
figure 3; convergence to the limiting incident-shock equilibrium conditions is apparent. The solu- 
tion is typical of the dissociative relaxation of a diatomic gas. The degree of dissociation is frozen 
(at a = 0) across the shock wave, and downstream of the shock all properties vary monotonically 
from values obtained assuming frozen composition to values corresponding to chemical equilibrium. 
Further, all the properties reach equilibrium on the same time scale. 

Variations in several end-wall properties as a function of time are shown in figure 4, for 
case 1. The most striking departure from the typical behavior exhibited in a steady incident- 
shock flow field is evident in the pressure history. Whereas the pressure behind the incident shock 
wave is nearly constant (fig. 3), the pressure behind the reflected shock wave first decreases in a 
comparatively short time span (~0.3 /xsec) and then increases over a much longer time period 
( — 1 00 jusec). The initial decrease may be attributed to the fast chemical relaxation occurring behind 
the reflected shock wave, whereas the subsequent increase is a consequence of the slower relaxation 
behind the incident shock wave (but upstream of the reflected shock wave). Note that, except for 
pressure, the equilibrium property values indicated in figure 4 are for reference only and do not 
represent asymptotic limits at the wall. 
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For case 1, the relaxation time behind the incident shock wave, r 2) is much longer than the 
relaxation time behind the reflected shock t 3 q (a graphic description of these characteristic times 
is shown in fig. 5). The gas in the reflected region thus quickly relaxes toward a state of local 
thermodynamic equilibrium corresponding to chemically frozen flow, with a=0, behind the 
incident shock and local chemical equilibrium behind the reflected shock. The pressure for such a 
condition is less than that corresponding to chemistry frozen at a = 0 behind both the incident and 
reflect shock waves since the effect of increasing density in the reflected region (due to dissociation) 
must be to decrease the reflected-shock-wave velocity and hence to decrease its corresponding 
pressure jump. 

The slow increase in pressure toward the reflected-shock “equilibrium” value is due to the 
nonequilibrium flow behind the incident shock wave. (A reflected-shock “equilibrium” value is 
defined here as the magnitude of a flow-field property when calculated assuming complete 
thermodynamic equilibrium behind both the incident and reflected shock waves; this is the value 
that applies in the flow far from the wall and is tabulated in tables of reflected-shock properties for 
real gases.) To a first approximation the pressure rise at the end wall occurs with the same time scale 
and in direct proportion to the density variation behind the incident shock wave. This finding is 
expected as the end wall effectively samples the incident-shock dynamic pressure, which is 
proportional to the density behind the incident shock wave. One consequence of this large pressure 
variation is that the end-wall pressure now becomes an attractive observable for reaction-rate 
studies. Similar behavior of the end-wall pressure was previously noted by Baganoff (ref. 10) and 
Johannesen et al. (ref. 1 1) for shock-wave reflection in vibrationally relaxing carbon dioxide. 

The temperature at the end wall has an initial value corresponding to a = 0 behind the 
reflected shock wave. As the relaxation proceeds, the temperature decreases rapidly toward the 
equilibrium value, and then increases slightly toward a final steady-state value. The presence of two 
time scales again is apparent as the temperature approaches a minimal value in a time that is much 
less than that required to approach the steady-state value. The time required for the temperature 
to approach its minimal value is related to the chemical relaxation time behind the reflected shock 
wave. This characteristic time is also evident in the relaxation of other variables, such as <x and S, 
which are sensitive to chemical processes. The density time history can be calculated from the 
results plotted in figure 4; its general character consists of a continuous increase in density with an 
overall time scale comparable to that of the slow increase in pressure. 

Another consequence of the nonsteady nature of shock-wave reflection in a relaxing gas is 
evident in the behavior of S, a, and T shown in figure 4 for long times on the end wall. All 
three of these properties approach steady-state values that are above their corresponding “equilib- 
rium” values. A plot of the density, on the other hand, would show a steady-state level below the 
equilibrium value. Although the pressure eventually must approach its equilibrium value, one should 
not expect the final steady-state solution for the other thermodynamic properties to correspond to 
the usual equilibrium values (since transport processes have been omitted here). Physically, the 
explanation for this result is that gas particles adjacent to the end wall are processed by a stronger 
reflected shock wave and have different time histories than those particles that are processed by the 
eventual equilibrium-strength reflected shock wave. Table 2 summarizes the end-wall values at 
41.4 /usee, at which time the solution for case 1 was terminated because of excessive computer 
time (approximately 6 hr on an IBM 7094 computer). A continuation of the computer solution 
would produce small additional increases in T, S, a, and p until the equilibrium value of pressure 
is attained. 
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TABLE 2.- COMPARISON OF END-WALL SOLUTIONS 


Solution 

p, 

atm 

T, 

°K 

5, 

cal/g— °K 

a 

pip i 

Nonequilibrium (41.4 psec) 

3.176 

3994 

2.558 

0.393 

63.80 

Equilibrium 

3.268 

3919 

2.501 

.342 

69.39 

Ratio = nonequilibrium/ 
equilibrium 

.97 

1.02 

1.02 

1.15 

.92 


To summarize the information presented in figure 4, the steady-state values on the end wall, 
for all of the thermodynamic properties except pressure, do not correspond to the usual equilib- 
rium values. Also, the relaxation of the flow appears to be controlled by two separate relaxation 
times. The first, t^c, is a consequence of the chemical relaxation behind the reflected shock 
wave. The second r 3 j, the time for the pressure to reach a steady-state condition, is controlled by 
the relaxation time behind the incident shock wave, t 2 ■ It is obvious that a steady-state solution 
cannot be attained behind the reflected-shock wave until the right-running characteristic that 
originates at the intersection of the reflected shock wave and the end of the incident-shock relaxa- 
tion zone, a loosely defined point, reaches the end wall. Thus r 3 j is approximately equal to r 2 . 

Flow-field properties for case 1 are presented in figure 6 as a function of distance from the 
end wall for discrete times after shock reflection. These plots are essentially “snapshots” of the 
various properties in the entire reflected-shock flow field taken at specified instants of time after 
shock reflection. Figure 6(a) shows the variation in pressure for times between 1 and 40 psec. 
Note that a plateau of near constant pressure develops within a few microseconds. The large pres- 
sure variation throughout the flow field for very short times, and near the shock for longer times, is 
associated with the chemical relaxation behind the reflected shock wave. The slow increase in the 
pressure level of the plateau can be attributed to the chemical relaxation behind the incident shock 
wave. The terminus of each curve in figure 6(a), away from the end wall, represents the pressure 
immediately behind the reflected shock wave at the specified instant of time. It can easily be shown 
that since a pressure gradient cannot be sustained in the flow field, the pressure must approach the 
equilibrium value throughout the flow field for long times after reflection. Thus, in a coordinate 
system that is fixed relative to the reflected shock, for long times after reflection the flow must 
show behavior similar to that of an incident-shock solution with upstream properties corresponding 
to the equilibrium conditions behind the actual incident shock wave. Downstream of the reflected 
shock wave, of course, there will ultimately be a steady nonequilibrium relaxation zone bounded by 
an over-growing region of gas in chemical equilibrium and at rest with respect to the shock tube. 

The spatial variation in entropy behind the reflected shock wave is shown in figure 6(b). The 
entropy at x = 0 quickly exceeds the equilibrium value (fig. 4), and figure 6(b) shows that it 
extends into the flow field in the form of an excess-entropy layer bounded by the end wall. As the 
shock moves away from the end wall, this layer becomes of nearly fixed extent at times greater 
than 1 0 psec and a steady-state entropy value is approached outside the layer that is asymptotic to 
the equilibrium value. The large entropy variations near the shock wave are a consequence of the 
chemical relaxation immediately behind the reflected shock. 
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Spatial distributions of the degree of dissociation and of temperature are shown in figure 6(c) 
and (d), respectively. Note that throughout the fully developed excess-entropy layer both a and T 
have values that are higher than their corresponding equilibrium levels. However, further from the 
end wall, steady-state plateaus begin to develop, at times greater than about 10 /zsec, wherein the 
properties approach their equilibrium values. Of course, large variations in a and T occur adjacent 
to the shock in the reflected-shock relaxation zone. As the flow properties upstream of the reflected 
shock become increasingly uniform, the relaxation profiles downstream of the reflected shock ap- 
proach steady-state distributions (as viewed in shock-fixed coordinates). After 40 Msec, for example, 
the reflected shock has propagated through the incident-shock relaxation zone a distance of about 
14 cm and the incident-shock properties have nearly reached their equilibrium values (fig. 3). 

Spatial variations in density are shown in figure 6(e) for the same time increments. The results 
indicate a general increase in density with time at all flow-field positions, although the final density 
throughout the entropy layer and near the shock front remains at a lower level than the equilib- 
rium values. Because the density is strongly coupled to pressure, no steady-state value is approached 
in the time span shown on the figure. This coupling between pressure and density has a simple 
physical interpretation. The. rising pressure compresses the gas, thus increasing its energy, but the 
heat-sink effect of the dissociation energy acts to hold the temperature and degree of dissociation 
nearly constant once “local” chemical equilibrium has been reached. Subsequent changes in pressure 
are thus accompanied by nearly proportionate changes in density. 

Measurements of density profiles have been used extensively to obtain rate information behind 
incident shock waves. The variation in density behind the reflected shock wave, at a fixed distance 
from the end wall, could also be used to infer reaction rates, but care must be exercised in inter- 
preting the data since the density profiles are different functions of time for different observation 
stations. The decrement in density associated with the excess-entropy layer should be detectable in 
a density measurement; in fact, Johannensen et al. (ref. 11) have obtained interferometer records 
exhibiting such a density decrease, partially masked by end-wall heat-transfer effects, for the 
vibrational relaxation of carbon dioxide behind a reflected shock wave. 

A boundary condition for “equilibrium” calculations behind reflected shock waves is that the 
flow velocity, relative to the shock-tube end wall, vanishes throughout the flow field. However, as 
mentioned earlier, for nonequilibrium calculations this boundary condition can be imposed only at 
the end wall. Spatial variations of velocity are shown in figure 6(f). Although a large gradient in flow 
velocity exists near the shock wave, except for small times the flow velocity approaches zero over 
a major portion of the flow field, as one would expect since the pressure gradient (spatial) is so 
small. 

Strehlow and Cohen (ref. 5) have shown that vibrational relaxation in nitrogen and oxygen 
causes an initial deceleration in the velocity of the reflected shock wave, relative to the shock tube. 
Similar behavior is shown in figure 7 for the case of dissociation of a diatomic molecule. The 
reflected shock velocity initially (at x = 0 ) has a value corresponding to a = 0 behind the 
reflected shock wave. However, the velocity quickly decelerates to a minimum value that is less 
than the equilibrium value and corresponds closely to the value computed under the assumptions of 
a. = 0 behind the incident shock wave and equilibrium chemistry behind the reflected shock wave. 
After reaching this minimum, the reflected-shock velocity increases slowly and should eventually 
reach the equilibrium value. 
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The foregoing figures present a complete description of the nonequilibrium reflected-shock 
flow field for case 1. Primary conclusions are: (1) two relaxation times control the flow behind the 
reflected shock; and (2) an excess-entropy layer is formed, near the end wall, in which the entropy 
and all other thermodynamic variables except pressure differ from their “equilibrium” values. The 
excess of entropy in this layer will not diminish, because within the assumptions of the present 
model there is no mechanism available to drive the entropy back to its equilibrium value. In the 
real case, of course, transport processes are present and would act to smooth out all existing flow- 
field nonuniformities. 

Capiaux and Washington (ref. 18) have investigated the analogous problem of nonequilibrium 
hypersonic flow of a dissociating gas past a semi-infinite wedge with a sharp leading edge. Their 
calculations also predict the formation of an entropy layer, located above the surface of the wedge, 
but they did not discuss the magnitude of the steady-state surface entropy in comparison to the 
equilibrium value. However, for their case, wherein the shock wave was convex, one also would 
expect that the surface entropy should exceed the equilibrium value appropriate for the deflection 
angle of the wedge. 


Cases 2 — 5 

With the general features of a nonequilibrium reflected-shock flow field understood, the 
results of cases 2 through 5 can be used in evaluating the influence of parameters such as initial 
pressure, reaction-rate constant, and shock velocity. We first consider the computed incident-shock 
flow fields for cases 2 through 5 (fig. 8) and compare them with the results of case 1 (fig. 3). 
Case 2 illustrates the effects of a fivefold increase in P y , with U and j3 held at their case 1 values. The 
primary effect is simply a decrease in the length of the relaxation zone, as might be expected 
because of the resultant increase in the magnitude of the rate equation (eqs. (7) and (1 1-13)). From 
a physical viewpoint, an increase in initial pressure increases the collision frequency of the molecular 
species, leading to a faster rate of reaction. Within the framework of the present model, some 
thermodynamic properties immediately behind the shock front (such as T, oc, p/p 1 , and PjP\) 
are independent of the initial pressure and hence take the same values for cases 1 and 2. The equilib- 
rium conditions are pressure dependent to some extent, however, because the equilibrium degree of 
dissociation is, a function of temperature and pressure. The net changes for case 2 compared to case 1 
equilibrium values are small decreases in a. and p/p 1 and a small increase in T\ the equilibrium 
pressure for case 2 is almost exactly five times that of case 1 . 

Case 3 demonstrates the effects of increasing the rate constant by a factor of 5 with the initial 
pressure and shock velocity held at their case 1. Since the rate parameter j3 does not enter the 
calculation of properties at any one flow-field point, but only affects the distance between adjacent 
flow-field points (eqs. (7)-(9), (17), and (18)), the flow-field profiles of case 3 are identical to 
those of case 1 except for a difference in the spatial scale of exactly a factor of 5. 

Cases 4 and 5 show the effects of changing the incident-shock speed U with and /3 fixed 
at their standard values (fig. 8). In case 4, as might be expected for the case of a stronger shock 
wave, the variations in all the properties are larger across the relaxation zone and the temperature 
level is generally higher. The increase in the equilibrium degree of dissociation is especially notable; 
this parameter is quite sensitive to the energy content of the gas, which in turn depends strongly 
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on the magnitude of U. A further point is that the overall relaxation length is seen to decrease 
with increasing shock speed, primarily because of the strong dependence of the rate constant on 
temperature. Case 5 represents the case of a weak shock wave, which causes very little dissociation. 
Accordingly, only small variations in property values are present in the incident-shock relaxation 
zone. The low temperatures are responsible for the great increase in the scale of the relaxation zone. 
Case 5 thus presents, effectively, the limiting situation of shock-wave reflection into a uniform 
nonreacting flow. 

Figure 9 compares end-wall properties plotted as a function of time after reflection for 
cases 2 through 5. With minor variations, discussed below, the results are similar to those of 
case 1 . Two characteristic times are present in each case, for example, and consistent deviations are 
found between the steady-state and equilibrium property values. To avoid excessive computer time, 
the solutions for all five cases were terminated prior to reaching a true steady state, but the final 
values of all properties can be estimated from the results shown. The computed pressure should 
approach asymptotically to the equilibrium value in all cases. 

The results for case 2, which illustrate the influence of increased initial pressure, differ from 
those of case 1 primarily in the time to reach steady conditions. The time scale of the entire 
solution is reduced because of the increased reaction-rate magnitude downstream of the reflected 
shock wave and the decreased spatial extent of the relaxation zone upstream of the reflected shock 
wave (fig. 8). For case 2, the general level of the pressure is higher than for case 1 , and because of its 
dependence on pressure, the entropy is smaller. The equilibrium property values for case 2 differ 
slightly from those for case 1 , but the deviations between steady-state and equilibrium values are 
nearly the same for both cases. 

The results for case 3 (fig. 9) are identical to those of case 1 (fig. 4) except for a factor of 
5 change in the time scale. Thus, the case 3 results confirm that the end-wall-properties solution, as 
well as the incident-shock flow-field solution, scales inversely with the magnitude of |3. Justifica- 
tion for this result follows from a consideration of the reflected-shock flow-field computation 
technique (appendix C) and an understanding of the incident-shock flow-field dependence on j3. 

The property-time histories of case 4 (fig. 9) show larger variations than all other cases, 
throughout the equilibration process, in all the properties plotted. This is the expected result for a 
stronger shock wave causing increased dissociation. In addition, the increase in shock strength leads 
to larger deviations between steady-state and equilibrium property values. The time to reach 
steady-state conditions is less than that for case 1 , as might be expected since the overall time is 
coupled to the relaxation time behind the incident shock wave r 2 , and r 2 is less for case 4 
than for case 1 . Another interesting feature of case 4 is the absence of an initial decay in end-wall 
pressure prior to the increase in pressure associated with chemical relaxation behind the incident 
shock. The initial period of near-constant pressure is actually a result of two competing effects. 
First, the chemical relaxation behind the reflected shock wave introduces a negative pressure gradient 
in time (with a time scale r 3 ^) as the flow relaxes toward chemical equilibrium. Second, the 
gradient in flow properties behind the incident shock wave induces a positive pressure gradient (with 
an approximate time scale r 2 ) on the end wall. Case 4 represents a solution wherein the initial 
values of these two gradients nearly cancel and no dip in pressure occurs on the end wall — that is, 
there is a strong interaction with the incident-shock flow field. A strong interaction is thus to be 
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expected whenever the relaxation times t 3 q and t 2 are not widely separate. In general, the dif- 
ference between r 3( o and r 2 lessens, and the extent of the interaction grows, with increasing 
shock strength. 

The results for case 5 (fig. 9) are typical of end-wall solutions for shock waves producing 
only a small degree of dissociation. In such cases, there is a weak interaction between incident- and 
reflected-shock flow fields, and the characteristic times t 3 q and t 3 j are completely distinct. 
The results shown in figure 9, however, are restricted to times of order t 3( j and do not approach 
the scale of r 2 or t 3 j- It appears that local chemical equilibrium is reached in approximately 
10 2 jusec, as evidenced by the behavior of all variables plotted, and the small increase in pressure 
shown for times greater than 10 2 psec indicates the onset of the influence of the incident-shock 
relaxation. The relaxation time r 2 is of the order of 10 4 psec for this case, so end-wall conditions 
at 10 2 psec should correspond closely to the limiting case of shock-wave reflection into a uniform, 
nonreacting flow. Note that an excess-entropy layer has already formed; the magnitudes of T,S , and 
a at 1 0 2 p sec are greater than their equilibrium values, so the formation of an entropy layer does 
not depend solely on the propagation of the reflected shock wave through a nonuniform upstream 
region. The effect of the incident-shock relaxation zone will be to greatly increase the spatial extent 
of the entropy layer, however, and to a lesser extent to increase the deviations between the steady- 
state and equilibrium wall-property values. 

It has been shown, in figures 4 and 9, that large variations in end-wall pressure may be present 
during shock -wave reflection. Pressure is a convenient quantity for measurement, but its suitability 
as a diagnostic tool in reaction-rate studies depends in part on the sensitivity of end-wall pressure 
to experimental parameters and the chemical reaction rate. For the purpose of further illustrating 
this sensitivity, the end-wall pressure histories for all five cases are replotted in figure 10, where, for 
convenience, the pressure coordinate has been normalized by the frozen (a = 0) reflected-shock 
pressure for each solution. From the results shown here (and in figs. 4 and 9), it is clear that the 
wall pressure is much more sensitive to the influence of chemical relaxation behind the incident 
shock wave than to the shorter-term effects of relaxation in the reflected region. Furthermore, it is 
obvious that the magnitude of the pressure variation following reflection increases with shock 
strength (compare the results for cases 4 and 5, for example) and that the time scale of the pressure 
history, for a fixed-speed shock, is sensitive to the parameters of initial pressure and rate constant 
(cf. cases 1-3). Thus it appears that measurements of wall pressure would provide an attractive 
method for studying reaction rates behind moderately strong or strong incident shock waves. 

Flow-field plots at discrete times after reflection are presented for cases 2 through 5 in 
figures 11 through 14, respectively. The format is the same as that used for case 1 in figure 4; 
since the results are generally similar in character to those of case 1, only new and significant 
features are discussed here. The distinguishing characteristic of the plots for cases 2 through 4 
is the reduction in the time required to reach near-steady conditions, the decrease resulting from 
the dependence of the relaxation time on the pressure, rate-constant multiplier, and shock velocity, 
respectively, for cases 2, 3, and 4. Of course, the reduction in relaxation time is also responsible 
for diminishing the entropy-layer thickness. One important result, seen from a comparison of 
figures 4 and 12, is that the reflected-shock flow-field distributions scale in a simple way with the 
value of j3. That is, all property values at the wall and at the shock front for case 3, at any time t, 
match exactly those for case 1 at the time 5 1 (recall that j3 = 1 for case 1 and j3 = 5 for case 3). 
At these times, the distances between the wall and the shock front for cases 1 and 3 differ by a 
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factor of 5, as do the gradients in each flow property. Thus, the entire shock-reflection solution, 
including the reflected-shock flow-field distributions, scales in a simple manner with the value of 
the rate-constant multiplier |3. 

Several interesting features appear in the results plotted for case 5 (fig. 14). Since the time 
scale of the plots ( t < 200 jasec) corresponds roughly to the relaxation time in the reflected region, 
the results for case 5 illustrate more clearly than those for the previous cases the progressive behavior 
of each property during this relaxation process. The results are typical of shock-wave reflection 
with chemical nonequilibrium downstream of the reflected shock and nonreacting flow upstream. 
The curves for pressure indicate that at about 100 /nsec a minimum is reached for the plateau pres- 
sure level, and for larger times an increase in pressure level is expected. The plots for the entropy 
distribution confirm that after a time of 200 /nsec the entropy layer is not yet fully developed, 
since no plateau in entropy is yet visible, and the entropy level remains higher than the equilibrium 
value. In fact, the entire reflected-shock flow field is still slightly out of equilibrium after 200 psec, 
as evidenced by the continuing changes in the degree of dissociation and temperature throughout 
the reflected region. After local equilibrium is achieved, the values of all properties will change 
further (although slowly) because of the effects of nonequilibrium behind the incident shock wave. 
These changes will be small because of the nearly negligible degree of dissociation throughout the 
incident-shock relaxation zone. 


Figure 1 5 compares results of cases 2 through 5 for the reflected shock velocity plotted as a 
function of distance from the end wall. The results are consistent with those of case 1 (fig. 7) with 
one exception: In case 5, the computations were stopped before a minimum in the shock-velocity 
history was reached. Further calculations should have shown a further decrease in reflected-shock 
velocity to a minimum value (less than the equilibrium value) followed by a slow, asymptotic rise 
toward the equilibrium value. 


CONCLUDING REMARKS 


In contrast to nonequilibrium incident-shock flow fields, large temporal variations in pressure 
may occur behind reflected shock waves. The magnitude of this pressure change and the sensitivity 
of its time scale to experimental variables make it a useful quantity for measurement in 
reaction-rate studies. 

The relaxation process behind a reflected shock wave exhibits two separate time scales. 
Properties that are sensitive to the chemical state, such as entropy, temperature, and the degree of 
dissociation, relax in a time that is controlled by the chemical reaction rate behind the reflected 
shock wave. Properties that are sensitive to fluid-dynamic phenomena, such as pressure and density, 
relax in a time r that is controlled by the time to obtain steady-state conditions behind the 
incident shock wave; t^j is usually much greater than T $c 
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An excess-entropy layer of gas is formed, adjacent to the end wall, in which all of the 
thermodynamic properties, except pressure, differ significantly from their usual equilibrium con- 
ditions as a steady-state solution is approached. 


Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, Calif., 94035, April 12, 1972 
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APPENDIX A 


SIMPLIFIED EXPRESSIONS FOR THE CHEMICAL EQUILIBRIUM 
CONSTANT AND THE CHEMICAL POTENTIALS 


The equilibrium constant for a dissociating, homonuclear diatomic gas may be written (ref. 1 6, 
p. 166) 


K _ [A ] 2 

Kc ' TaT) 


eq 


(q^ivy 

N{QA 2/E) 


moles/ cm 3 


(Al) 


where A and A 2 denote atomic and molecular species, respectively, Q ^ and Q^ 2 are the atomic 
and molecular partition functions, V is the volume of the system and N is Avogadro’s number. 


For the atoms, the complete partition function is simply the product of the translational and 
electronic partition functions 

C 4 = Q fr < 4 e ^ dl2T (A2) 

The exponential factor results from choosing as the reference energy level the ground vibrational 
and electronic state of the molecule. The translational partition function is given by 

qA = ^2n m^ ry 2 v (A3) 


where m is the mass of the atom, and k and h are the Boltzmann and Planck constants, 
respectively. The electronic partition function is given by 



(A4) 


where the g\ are degeneracy factors for the levels of energy characterized by the temperatures 
®el,/ • 

For the molecules, the complete partition function is the product of the translational, 
rotational, vibrational, and electronic partition functions 


qA 2 = qA 2 Q A 2 qA. 2 A A 2 

^ ^tr ^rot ^vib el 


(A5) 
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The translational and electronic partition functions take the same form as expressed above for the 
atoms, but are evaluated with the molecular mass (2m) and molecular electronic-level parameters. 
The rotational and vibrational partition functions are (ref. 1 6) 



(A6) 


and 



(A7) 


where d r and 6 V are characteristic rotational and vibrational temperatures. 

A simplified expression for the equilibrium constant in oxygen can now be formed by sub- 
stituting appropriate values for the constants involved in the expressions above. Assuming approxi- 
mate values for the electronic partition functions of 8.8 and eS? = 3.3 (reasonably valid 
over the temperature range of interest in this work), one finds that 

K C (T) = 0.666 (T) l/2 (l - e~ dvlT ^j e~ 6d ^ T moles/cm 3 (A8) 

The chemical potential of each component in a mixture of thermally perfect gases is given 
by (ref. 1 6) 


m = 



cal/g 


(A9) 


where Qi is the complete partition function of species i, Wj is the molecular weight, and N{ is 
the number of particles of species i in the system. 


From the above discussion of partition functions, and again choosing the ground molecular 
state as the reference energy level, one can immediately write 



e ~ed/2T 


(A10) 


and 


qAi _ ^nkT_ 
A U 2 ” h 2 



(All) 


Simplified expressions for oxygen are again formed by substituting values for the constants involved 
and noting that F/A^ = kT/p^ (where is the partial pressure of the atoms, in dynes/cm 2 ): 


Q° 
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(A 12) 



and 
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By expressing the partial pressures in units of atmospheres (rather than dynes/cm 2 ), one obtains 
final simple expressions 


Ho = - (*7716) In [l4.419(r s/2 lp 0 ) e 6dl2T ] cal/g (A14) 

and 

Ho 2 = - (*7732) In [3.714(T 7 / 2 /p 02 ) (l ~ e~ 6vlT ) ] cal/g (A15) 

The partial pressures are related to the total pressure P by 

p °‘rh F <A16) 


and 
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1 - a 
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APPENDIX B 

METHOD OF CHARACTERISTICS FOR ONE-DIMENSIONAL 
MOTION OF A DISSOCIATING, HOMONUCLEAR DIATOMIC GAS 


The characteristic equations for pressure and mass velocity in a reacting gas are formulated in 
the same manner as that used for the case of a nonreacting gas — by manipulating general thermo- 
dynamics relations together with the conservation equations. For a homonuclear diatomic gas 
undergoing the reaction 


kd 

A 2 + X ^ 2A + X 
k r 


where X is any collision partner, the pertinent relations (see earlier text for notation) are as 
follows: 


Conservation of mass 


Conservation of momentum 


Equation of state 




= 0 



+ = ° 

ox 


p = p(P,S,a) 


Rate of entropy production owing to chemical nonequilibrium 


(Bl) 


(B2) 


(B3) 


DS 1 v~' 1 / \ 
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i 

Rate of production of atom mass density owing to chemical reaction 


w 4 = p— = W \kdlA 2 ][X] - k r [A] 2 [X]\ 


(B4) 
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Speed of sound at constant chemistry 


= (<*) 

\9p / S,a 


(B6) 


Equation (B3) may be written in differential form as 


Dp / 3p\ DP + /dp\ DS /9p\ 

Dt \dP )s a Dt \9 S/u,pDt \da)pgDt 

However, from the combined first and second laws of thermodynamics for this reacting system 
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In addition, h may be written in terms of P , a , and p so that 

* = (^L ^ + (^L * + (^k 


(B9) 


(BIO) 


(Bll) 


An inspection of this relation and the preceding differential expression for h allows one to write 

(PA~PA l 2 ) “ (9/2/9«)p,p 


( 


da ) P,S 


(dh/dp)p a 


(B12) 


Equation (B7) may now be rewritten as 


Dp = (*p\ 21 + t /Jp\ RS (PA^PAt) ~ (dhfdoQP'P 
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Upon substituting equations (B4), (B5), and (B6), one obtains 


Dp _ _1 DP _ (bh/ba)p rP wa 
D t a 2 dt (bh/bp)p >0L p 


(B14) 


The characteristics equations can now be formed by substituting equation (B 1 4) into equation 
(Bl), multiplying equation (B2) by 1/a, and first adding and then subtracting the resulting forms of 
these two equations. The results are 


and 
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These two equations determine the rate of change of pressure and mass velocity along the right- 
and left-running characteristic lines specified by dx/dt = u ±a. 


The characteristics equations can be cast into an abbreviated form with the notation 
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so that 


and 


1 

„2 8t a bt p 


(B20) 


1 8~P _ l_b~u_ _ G_ 
p a 2 bt a bt p 


(B21) 


24 



Thus 5 + /Sf and 8~/8t correspond to differentials along the right- and left-running characteristics, 
respectively. Eauations (B20) and (B21) can be used to calculate changes in P and u throughout 
characteristic network, once an appropriate form for the term G is obtained. Note that G = 0 
when = 0 in which circumstances the forms of equations (B20) and (B21) reduce to the 
standard equations used for isentropic flow. In the case of chemical equilibrium, the sound speed 
a takes on its equilibrium value. 

An explicit form for G is established by utilizing the thermal and caloric equations of state 
(eqs. (5) and (6) in text). From these expressions, one can write 


Y) .&) (f) -v(-) 

^ P/P,a \^ T )p,ol \3p /P,a Pi \ p ) (B22) 


where 


pf 



is the specific heat at constant pressure and composition; that is, 
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One can also write 


cal/initial mole - °K (B23) 


where 
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Substituting the above results in equation (B17), one finds that G is simply given by 


G 



R o»/r\ c p f - j w A 

e e v /T_ |/ 1 +aj Cpf 


where Cpf is specified by equation (B23). 
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APPENDIX C 

NETWORK CONSTRUCTION TECHNIQUE 


As in all characteristic solutions, the 
present solution proceeds from a known 
set of points in the flow field. Disregarding 
the starting approximation for the time 
being (appendix D), the solution proceeds 
from the known set of properties along the 
left-running characteristic ray 1-3, as shown 
in sketch (a). The properties on a subsequent 
ray, ray 4-6, are obtained from those on 
the initial ray by calculating three types of 
points: (1) end wall, (2) interior, and (3) 
shock wave. 



Sketch (a).— Characteristic network construction. 


End-Wall Point 

The procedure required to obtain values of the flow-field variables at an end-wall point, 
say point 4, is as follows: 

1 . The location of point 4 in the xt plane is established geometrically by the intersection 
of the right-running characteristic from point 2, slope dx/dt = (u + a) 2 and the end wall at x = 0. 

2. The pressure P 4 is obtained by solving the finite-difference form of equation (23a): 

_±\ P* ~ P 2 _ 1 *2 = 

pa 2 J 2 t 4 - h a 2 t 4 - t 2 

since the velocity at the end wall u 4 is zero. 

3. The values for the degree of dissociation a 4 and the entropy S 4 are found by solving 
the finite-difference forms of equations (29) and (30) along the end wall, which is the location of a 
particular particle path. To minimize accumulative errors, we first divide the overall time interval 
t 4 — ti into 20 equal subintervals (throughout which the pressure is assumed to vary linearly 
from Pj to P 4 ), and then apply equations (29) and (30) to each subinterval. The appropriate 
forms of equations (29) and (30), respectively, are 



<H ~ «f-i 

H ~ H-\ 



(C2) 
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and 


Sj - -Si-i r 1 


(C3) 


where i - 1 represents the point along the particle path at which the calculations have been 
completed and the properties are known, and i is the subsequent point at which time and pressure 
are known and all other properties are to be determined. 

4. Using the values of P, a, and S at each point i, the temperature can be found by satisfying 
the following equation for entropy (as a function of T, a, and P (in atm)) implicitly: 


c _ m +a) 

b W 


E^br) " E N k lnN k -InP 

. k V /k k 


(C4) 


where N k are mole fractions of the species involved, 

2a 


No = 


1 + a 


and 


No, = 


1 - a 
1 + a 


and ( S°/R) k is the dimensionless entropy of species k evaluated for a pressure of 1 atm; that is, 

In ^14.419 T 5 l 2 e~® d ^ 2T ^j 


(—) 

\R /, 


= 2 

2 + 2T + 


and 


/ S°\ _ 7 

0^> T , 

-t- § 

~ 3.7 14 T 7/2 ~ 

\ R )o 2 2 

r(i-^) 

j - ‘~ evlT . 


The calculations continue until a, S, and T are found at point 4. Other desired thermodynamic 
quantities can now be computed directly in terms of the known properties. 


Interior Point 

The solution of an interior point, say point 5 in sketch 1 , proceeds in much the same manner 
as described for an end-wall point. 

1. Point 5 is located geometrically by the intersection of a right-running characteristic 
from point 3, slope dx/dt = (u + «) 3 . and the left-running characteristic from point 4, slope 
dx/dt = ( u ~a) 4 . 

2. The pressure and velocity at point 5, P 5 and u 5 , are found by solving the following two 
equations simultaneously: 
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and 



P 5 — P 3 + 1 u 5 — u 3 _ /G \ 
t s - h a 3 t s - t 3 \p l 



1 u 5 - U 4 

a 4 t 5 ~ ^4 



(C5) 


(C6) 


3. The degree of dissociation and the entropy at point 5 are found by the technique described 

in step 3 for the end-wall point — that is, by applying the finite-difference forms of equations (29) 
and (30) to 20 subintervals along the particle path passing through point 5. The origin of this 
particle path, point a, is found by searching along the characteristic segment between points 2 and 3 
(using linear interpolation) for a value of particle velocity that when averaged with the known 
velocity at point 5, corresponds to passage through point 5. The values of P , T , and a at point a 
are found by linear interpolation between the known properties at points 2 and 3. Values of other 
necessary quantities at point a, such as p, S, w Q , po> and po 2 > are then computed from the 
interpolated values of P, P, and a . 


4. The temperature at each subinterval point along the particle path is found from a 
knowledge of P, a, and S at the corresponding point, as described for the end-wall point. 


Shock-Wave Point 

Values of the flow-field variables at a shock-wave point are obtained as follows: 

1. Coordinates of a shock-wave point, point 6, are found from the intersection of a left- 
running characteristic from point 5, slope dxjdt = ( u~a ) 5 , and the reflected shock wave from 
point 3. The velocity of the reflected shock wave is assumed to be constant between points along the 
shock wave trajectory, and is allowed to change discontinuously at each shock-wave point. 

2. Flow properties in front of the reflected shock wave are found from the appropriate 
position in the incident shock-wave solution — that is, at a distance s behind the incident shock 
where (for point 6) 


$ 6 = -*6 + Ut 6 

3. Pressure and particle velocity immediately downstream of the shock, at point 6, are 
found by solving 

/ 1 \ P 6 - P 5 1 u 6 — u 5 _ /g \ (C7) 

\pa 2 3 J s t 6 - t 5 a 5 t 6 - t 5 \p J 5 

simultaneously with the Rankine-Hugoniot jump conditions across the reflected shock wave. The 
only boundary condition imposed on the solution of this point is that the degree of dissociation 
remains invariant across the reflected shock wave. The additional thermodynamic properties needed 
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and a new value for the reflected shock speed follow directly from the conservation equations 
(applied across the reflected shock) and the state relations. 


Methods of Minimizing Error 

The flow-field solution is continued by successive computation of rays. Although the 
description herein is for a fixed number of points (each subsequent ray increases the number of 
points per ray by one), additional points may be added (by interpolation) along each successive ray 
to keep the mesh sufficiently small to preserve accuracy. Furthermore, to obtain good accuracy and 
stability in the region immediately behind the reflect shock wave, where very large property 
gradients are present, a fine mesh is incorporated into this region of the network (see appendix E). 
The number of points in this fine mesh is variable and depends on the magnitude of the property 
gradients. 

In addition to the procedures described above for calculating the three types of points, the 
solution at each point in the entire flow field is iterated once. That is, in order to minimize errors 
due to large property gradients, it is more efficient to iterate the solution of each new flow-field 
point once than to use a finer network of points and not to iterate. The technique simply involves 
averaging the properties of a newly calculated point with those of the known or reference points, 
assigning the average values to the coordinates of the known points, and recomputing the position 
and properties of the new point. 

Iteration at each flow-field point and the use of a fine mesh scheme near the shock wave have 
proven effective methods of minimizing errors (e.g., the so-called entropy diffusion effect) which 
accumulate or propagate in the solution. 
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APPENDIX D 


STARTING APPROXIMATION 


To obtain the flow-field properties along the initial left-running characteristic, ray 1-3 in 
sketch (a), properties are first calculated at the end-wall points 0 and 1. Properties at point 0, 
shown as the origin of the xt diagram, are obtained from the standard shock-jump relations for 
shock-wave reflection at a plane, nonconducting wall. As is appropriate for the instant of reflection, 
a = 0 behind the reflected shock wave. Properties at the second point, point 1 , can be approximated 
in a number of ways. However, it is desirable to have a good approximation so as to minimize the 
perturbation introduced into the resulting solution. Since the velocity on the end wall vanishes and 
heat-transfer effects are ignored, it is reasonable to assume that the enthalpy along the end wall is a 
constant for a short time — that is, h x = h 0 . The degree of dissociation at point 1 can then be 
found in terms of the temperature from 


where 


= 


ho - (7/2) - g Vl 

(3/2) RT x - e Vl + D 


e Vl = 


Rdy 
fivlTi _ 


(Dl) 


(D2) 


is the vibrational energy at point 1 . The temperature has been used as the independent variable, 
primarily for convenience, and is evaluated at point 1 using 


Ti = T 0 - AT 


(D3) 


where AT is a specified initial perturbation. The third property required to define the nonequilib- 
rium state at point 1 is obtained by assuming, consistent with the constant-enthalpy approximation, 
that P x = P 0 . 

The time duration of the starting approximation At , defined as 

2a!! 

At = 1 1 “ to (doi/dt ) i + (da/dT ) 0 (D4) 

is, of course, required to be very small in comparison with the total time scale of the relaxation 
process. Thus, the assumptions of constant enthalpy and pressure during this interval should not 
seriously affect the resultant solution. The starting-time interval is controlled by the magnitude of 
the initial temperature perturbation AT . The magnitude of AT for each case was selected to 
yield an optimum mesh spacing, by a trial-and-error procedure, for the resulting solution. Although 
the initial gradients of flow-field variables are somewhat insensitive to AT, the solution can become 
unstable at large times if AT is chosen too large. On the other hand, if AT is too small, excessive 
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computer time will be required to obtain a solution. A compromise for the cases considered here 
put AT in the range 1-0° < AT < 400° K. 

The remaining points shown along ray 1-3 are obtained by first calculating the properties at 
point 3 using the procedure for shock-wave points described in appendix C. The position of point 3 
is located by the intersection of the reflected shock wave from point 0, slope dxfdt = ~V 0 and 
the left-running characteristic from point 1 , slope dx/dt = — a t . Linear interpolation between points 
1 and 3 is used to obtain the properties at the intermediate point 2. The ray 1-3 is now ready for 
use as a known ray for the subsequent flow-field calculations. 
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APPENDIX E 


FINE-MESH SCHEME AND ROLLBACK TECHNIQUE 

To maintain sufficient accuracy and stability in the region immediately behind the reflected 
shock wave, where very large gradients are present in the thermodynamic properties, a fine mesh is 
incorporated into this region of the network. The number of points within this fine mesh is 
dependent on the magnitude of the temperature gradient adjacent to the reflected shock wave. 

Sketch (b) is a diagram of the points 
involved in a typical, narrow segment of 
the fine-mesh network. For clarity, left- 
running characteristics calculated prior to 
ray 1-6 are not shown; this discussion 
concerns the logic of the fine-mesh net- 
work that follows after points 1-6 are 
known. The points are labeled in the order 
that they are calculated. 

The calculation technique for obtain- 
ing the position and properties of points 7 
through 1 0 is described in appendix C. The 
next step is to select the number of points 
to be added to the basic network. In the 
present study, one new point is added 
along the line segment between the shock- 
wave point and the preceding point (be- Sketch (b). -Fine-mesh characteristic network, 

tween points 6 and 5) for each 400° K of 
temperature difference. For example, if 

the temperature difference between points 6 and 5 is less than 400° K, a single new point is added, 
as shown by point 1 1 in sketch (b). If the temperature difference is between 400° and 800° K, 
then two new points (equally spaced) are added; for example, see points 25 and 26 along the 
next continuous left-running characteristic in the sketch. 

The properties P 9 u , T, and a at point 1 1 are found by interpolating between the known 
property values at points 4, 5, and 6 (using a quadratic fit for the property values as a function of 
distance along the characteristic). Other needed properties are calculated directly from these inter- 
polated property values. The number of points added (always equally spaced) along the ray from 
point 5 to 10 is the same as added between points 5 and 6 — in this case, one. The properties P 
and u at that point, labeled 12, are established by linear interpolation. Other properties are com- 
puted with the standard interior-point technique involving streamline integration and linear inter- 
polation between points 4 and 5. 

On the basis of the known properties at points 5, 6, 11, and 12, the location and properties 
at points 13 and 14 follow directly from the standard interior-point routine. Point 15 is then 
established using the conventional shock-point routine. The auxiliary ray 12-15 now enables 
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calculation of points 1 6 through 1 9 (using standard routines) with an improved accuracy over that 
attainable without the auxiliary ray. 

The calculation of ray 20-43 is intended to represent the case where the temperature difference 
between points 1 9 and 1 8 is between 400° and 800° K. The order of the calculations is self explan- 
atory. The techniques are the same as described above for the case of a single added point adjacent 
to the shock wave. In this case, however, an additional point, point 44, is added midway between 
points 4 1 and 42 (by interpolation) to maintain more uniform mesh spacing in the next calculated 
characteristic. In fact, to maintain uniform spacing, only points 21-24, 40, 44, and 43 are used 
during calculation of the subsequent complete left-running ray. 

It should be noted that, in addition to the calculations described here, each of the points in the 
fine-mesh part of the network is iterated once using the technique discussed in appendix C. 

Since each new continuous left-running characteristic involves one more computing point than 
the previous ray, it is obvious that the number of computing points can rapidly become unmanage- 
able from the viewpoints of computer storage, time to achieve a solution, and printout-statement 
volume. To keep the number of calculations to a reasonable level, a rollback technique is employed. 
In short, each time a ray reaches 175 points, the ray is redivided into 1 50 equally spaced points with 
the properties at each new point determined by linear interpolation. 
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IDEAL FLOW NONEQUILIBRIUM FLOW 





DISTANCE FROM INCIDENT SHOCK, Cm 


Figure 1 Comparison of reflected-shock wave x-t Figure 3.— Incident-shock flow-field profiles for case 
diagram for ideal and nonequilibrium flow. 



Figure 2.— Characteristic network construction. 
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Figure 4.— End-wall properties as a function of time 
after reflection, case 1 . 




Figure 5 Graphic description of controlling times for 
reflected-shock-wave relaxation. 



DEGREE OF DISSOCIATION , a 



(a) Pressure. 



(d) Temperature. 
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DISTANCE FROM END WALL , cm 


(b) Entropy. 
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(c) Degree of dissociation. 



DISTANCE FROM END WALL , cm 


(e) Density. 



(f) Flow velocity relative to end wall. 


Figure 6.— Flow-field properties as a function of distance from the end wall for discrete times 

after reflection, case 1 
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DISTANCE FROM END WALL, Cm 


Figure 7.— Reflected-shock velocity relative to 
end wall, case 1. 



(a) Case 2. 


(c) Case 4. 



(b) Case 3. 


(d) Case 5. 


Figure 8.— Comparison of incident-shock flow Fields, cases 2-5. 


39 




(a) Case 2. 


(c) Case 4. 



(b) Case 3. 



(d) Case 5. 


Figure 9.— Comparison of end-wall properties as a function of time after reflection, cases 2-5. 



TIME AFTER REFLECTION, /XSeC 


Figure 10.— Comparison of end-wall pressure 
histories, cases 1-5. 
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Figure 1 1 Flow-Field properties as a function of distance from the end wall for discrete times 

after reflection, case 2. 
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DEGREE OF DISSOCIATION , a 
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Figure 12.- Flow-field properties as a function of distance from the end wall for discrete times 

after reflection, case 3. 
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DISTANCE FROM END WALL , Cm DISTANCE FROM END WALL , Cm 

(a) Case 2. (c) Case 4. 



(b) Case 3. 


(d) Case 5. 


Figure 15.— Comparison of reflected-shock velocities, cases 2-5. 
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